/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// The main program for the radiative transfer code. \ingroup mcrx

/** \defgroup mcrx Radiative-transfer calculation (mcrx)

    The mcrx executable is the actual radiative-transfer program.  It
    does the ray tracing, both for the non-scattering stage and the
    scattering stage, and then the necessary postprocessing:
    calculating attenuation, interpolating over wavelength, adding the
    IR template, etc. 

    The configuration options and output file format for mcrx are
    described on the \ref mcrx-readme page.
*/

// $Id$

#include "config.h"
#include "mcrx.h"
#include "mcrx-stage.h"
#include <string>
#include "CCfits/CCfits"
#include <memory> 
#include "single_lambda_emergence.h"
#include "emergence-fits.h"
#include <sys/resource.h>
#include "boost/thread/thread.hpp"
#include "hpm.h"
#include "misc.h"
#include "fcntl.h"
#include "threadlocal.h"
#include "mcrxversion.h"

#ifdef WITH_CUDA
#include "cuda_grain_temp.h"
#endif

#ifdef WITH_AREPO
#include "mcrx-arepo.h"
#endif

namespace mcrx {
  template<typename> class adaptive_grid;
  template<typename> class arepo_grid;
};

using namespace CCfits;
using namespace std;


/// \name SCATTERING_LAMBDAS HDU column names
///@{
/// Entry in lambda table.
const string mcrx::Mcrx::sllec = "entry";
/// Requested wavelength.
const string mcrx::Mcrx::slrlc ("rlambda");
/// Normalization.
const string mcrx::Mcrx::slnc ("normalization");
/// Number of rays.
const string mcrx::Mcrx::slnrc ("nrays");
///@}

namespace mcrx {
  /// Handy typedef for a pair of floats. \ingroup mcrx
  typedef std::pair<mcrx::T_float, mcrx::T_float> float_pair;
};
namespace std {
  /** Input operator for the float_pair, used for using istream_iterators
      when we read the camera positions file. \ingroup mcrx */
  std::istream& operator>> (std::istream& s, 
			    std::pair<mcrx::T_float, mcrx::T_float>& p)
  {
    s >> p.first;
    if (s) {
      s >> p.second;
      if (!s)
	cerr << "operator>> (istream&, float_pair&): Could not read both components!"  << endl;
    }
    
    return s;
  }

  /** Element-wise multiplication operator for the float_pair, used
      for transforming the input when we read the camera positions
      file. \ingroup mcrx */
  std::pair<mcrx::T_float, mcrx::T_float>
  operator*(const std::pair<mcrx::T_float, mcrx::T_float> l, 
	    const std::pair<mcrx::T_float, mcrx::T_float> r)
  {
    return std::pair<mcrx::T_float, mcrx::T_float>(l.first*r.first, 
						   l.second*r.second);
  }
};


/** Initializes the mcrx output file.  This function copies the
    necessary info from the input file to the output file and creates
    extensions describing the radiative transfer wavelengths and the
    emergence cameras.  */
void mcrx::Mcrx::initialize_output_file ()
{
  if (terminator ())
    return;

  cout << "\n  *** ENTERING INITIALIZATION STAGE ***\n" << endl;

  // First see if the output file exists
  const string output_file_name = 
    word_expand (p.getValue("output_file", string ())) [0];
  auto_ptr<FITS> output;
  try {
    cout << "Trying to open output file " << output_file_name << endl;
    output.reset(new FITS (output_file_name, Read));
    cout << "Checking if initialization stage already complete " << endl;
    // also check if the keyword mcrx_init is defined, in that case
    // this stage is complete and we can just return
    bool init;
    output->pHDU ().readKey("MCRX_INIT", init );
    if (!init) {
      cout << "init false" << endl;
      throw 0; // just throw whatever
    }
  }
  catch (...) {
    // the output file either didn't exist, or the init stage was not
    // completed.  In this case we simply re-create the file.
    cout << "Initializing output file " << output_file_name << endl;

    cout << "opening input file "  << '\''
	 << word_expand (p.getValue("input_file", string ())) [0]
	 << '\'' << endl;
    auto_ptr<FITS> input 
      (new FITS (word_expand (p.getValue("input_file", string ())) [0], 
		 Read));

    cout << "copying primary header" << endl;

    // "copy primary HDU" constructor
    output.reset(new FITS ("!" + output_file_name, *input));
    cout << "writing primary keywords" << endl;
    // Write the primary keywords
    output->pHDU().addKey("FILETYPE", string ("TRANSFER"),
			  "File contains output from a radiative transfer calculation");
    output->pHDU().addKey("DATATYPE", string ("GRID"),
			  "Data is on a grid");
    output->pHDU().addKey("TRANSFERCODE", string ("MCRX"),
			  "See HDU MCRX for radiative transfer setup");

    
    {
      // Create scattering_lambdas HDU
      cout << "Creating SCATTERING_LAMBDAS HDU" << endl;
      Table* scattering =
	output->addTable("SCATTERING_LAMBDAS", 0,
			 vector<string> (), vector<string> (), 
			 vector<string> (),
			 BinaryTbl); 
      scattering->writeComment("This HDU contains the number of rays shot at different wavelengths.  The first entry, with lambda = -1 is the no scatter run.");
      scattering->writeComment("Column " +slrlc+ " is the requested wavelengths in the configuration file.  These may differ from the actual wavelengths used because we use the closest entry in the LAMBDA table");
      // we need to get the units here
      scattering->addColumn(Tint, sllec, 1, "Entry in LAMBDA table");
      ExtHDU& input_lambda_HDU = open_HDU (*input, "LAMBDA");    
      scattering->addColumn(Tdouble, slrlc, 1,
			    input_lambda_HDU.column("lambda" ).unit(), 8);
      scattering->addColumn(Tlogical, "line", 1);
      scattering->addColumn(Tdouble, slnc, 1, "Same as corresponding image values", 8 );
      scattering->addColumn(Tlong, slnrc, 1, "" );
      scattering->column(slnrc).addNullValue (long (- 1));
      scattering->addKey("poly",true,"This is a polychromatic simulation");
    }
    
    // we want to copy all information from the input grid file EXCEPT
    // GRIDDATA, PARTICLEDATA and INTEGRATED_QUANTITIES.
    cout << "Copying HDU's from input file" << endl;
    try {
      
      const CCfits::ExtMap extensions = const_cast<const FITS*>(input.get())->extension();  
      for (multimap< string, ExtHDU*>::const_iterator i = extensions.begin();
	   i != extensions.end(); ++i) {
	const string ext_name (i->first); 
	if ((ext_name.find(string ("GRIDDATA")) != string::npos) ||
	    (ext_name.find(string ("PARTICLEDATA")) != string::npos) ||
	    (ext_name.find(string ("INTEGRATED_QUANTITIES")) != string::npos))
	  continue;
	cout << ext_name << endl;
	
	output->copy(input->extension(ext_name) );
      }
    }
    catch(...) {
      cerr << "mcrx::initialize_output_file: Error copying HDU's, output file is incomplete!"  << endl;
    }
   
    {
      // Generate wavelength grid for the intensity stage.
      // We use a logarithmic wavelength grid covering the same range as
      // the original one
      array_1 lambda;
      ExtHDU& lambda_hdu = open_HDU (*output, "LAMBDA");
      read(lambda_hdu.column("lambda"),lambda);
      const int n_lambda = lambda.size();

      const int n_lambda_int =
	p.getValue("n_wavelengths_intensity", int());
      if(n_lambda_int>n_lambda) {
	// if this is true, things will crash because it reads lambda
	// as being the length of the table and if n_lambda_int is
	// greater then the length of the table is set by that and the
	// lambda vector will contain junk at the end. We would have
	// to add a nlambda keyword and only read in that many (like
	// we do for the lambda_intensity) but it's a waste because it
	// doesn't make sense to have more wavelengths for intensity
	// anyway so we just exit.
	cerr << "You have requested higher wavelength resolution in intensity calculation.\nThis will actually make the code crash and is almost certainly not what you want, so\nI'm going to stop here."<<endl;
	exit(1);
      }

      const T_float dr=
	log10(lambda(lambda.size()-1))-log10(lambda(0));

      array_1 lambda_int(n_lambda_int);
      lambda_int = 
	pow(10, blitz::tensor::i*dr/(n_lambda_int-1)+log10(lambda(0)));

      // make sure the end points are *identical* regardless of truncation
      lambda_int(0) = lambda(0);
      lambda_int(n_lambda_int-1) = lambda(lambda.size()-1);
      // write to LAMBDA HDU
      lambda_hdu.addColumn(Tdouble, "lambda_intensity", 1,
			   lambda_hdu.column("lambda" ).unit());
      lambda_hdu.addKey("NLAMBDA_INTENSITY",n_lambda_int,
			"Number of wavelengths in intensity grid");
      make_onedva(lambda_int).write_column(lambda_hdu, 
					   lambda_hdu.column("lambda_intensity"),
					   1,0, n_lambda_int);
    }

    // set up units and add them to the preferences so they get
    // written to the MCRX HDU
    ExtHDU& pd_hdu = open_HDU(*input, "PARTICLEDATA");
    units["length"] = pd_hdu.column("position").unit();
    units["time"] = pd_hdu.column("age").unit();
    units["mass"] = pd_hdu.column("mass").unit();
    for(T_unit_map::const_iterator i=units.begin(); i!=units.end(); ++i)
      p.setValue(i->first+"unit", i->second);

    cout << "Writing configuration parameters" << endl;
    // write mcrx parameters to HDU "mcrx"
    write_parameters (*output);

    // Setup emergence
    setup_emergence (*output);

    // get seed and initialize the random number state
    cout << "Initializing random number generators" << endl;
    const int seed = p.getValue("seed", int ());
    mcrx::seed (seed);
    const int n_threads = p.defined("n_threads")?
      p.getValue("n_threads", int ()): 1;
    rng_states = generate_random_states (n_threads);
    write_rng_states (*output);
  
    // and when we are completely done, set the "stage complete" keyword
    output->pHDU().addKey("MCRX_INIT", true, "Initialization stage complete");
  }//catch

  cout << "Output file initialization complete" << endl;
}


/** Writes the configuration parameters to the MCRX HDU.  */ 
void mcrx::Mcrx::write_parameters (CCfits::FITS& file) const
{
  Table* parameters = file.addTable("MCRX", 0);
  p.write_fits(*parameters);
  write_history (*parameters);
}
  

/** Writes the CAMERAi-PARAMETERS HDU's to the output file.  The only
    way to do this is to actually create the emergence object which
    then writes itself.  This is the only time the camera
    configuration parameters are used, in the actual stages the
    emergence object is loaded from the information in the output
    file.  */
void mcrx::Mcrx::setup_emergence (CCfits::FITS& file) 
{
  // the "CAMERAi-PARAMETERS" are common to both nonscatter add
  // scatter stages, so we write them here.  The only way to do this
  // is to create an emergence object which seems kind of wasteful,
  // but whatever...
  cout << "Writing emergence CAMERAi-PARAMETERS HDU's" << endl;

  const std::string camtype = p.defined("camera_projection") ?
    p.getValue("camera_projection", string ()) :    
    rectilinear_projection::typestring_;

  // See if we are using Arepo and should translate the origin
  ExtHDU& sfrhist_hdu = open_HDU(file, "SFRHIST");
  string nbodycod;
  sfrhist_hdu.readKey("NBODYCOD", nbodycod);
  vec3d translate_origin(0,0,0);
  if(nbodycod=="arepo") {
    try {
      sfrhist_hdu.readKey("TRANSLATE_ORIGINX", translate_origin[0]);
      sfrhist_hdu.readKey("TRANSLATE_ORIGINY", translate_origin[1]);
      sfrhist_hdu.readKey("TRANSLATE_ORIGINZ", translate_origin[2]);
    }
    catch (HDU::NoSuchKeyword&) {
      translate_origin=0;
    }
    cout << "Translating camera origin to " << translate_origin << endl;
  }

  auto_ptr<full_sed_emergence> e;
  const int npixels = p.getValue("npixels", int ());
  if (p.defined("ntheta") && p.defined("nphi" )) {
    const T_float fov = p.getValue("camerafov", double ());
    const int ntheta = p.getValue("ntheta", int ());
    const int nphi= p.getValue("nphi", int ());
    const bool exclude_south = p.defined("exclude_south_pole") &&
      p.getValue("exclude_south_pole", bool ());
    const T_float distance = p.getValue("cameradistance", double ());
    e.reset(new full_sed_emergence (ntheta, nphi, distance,
				    fov, npixels, exclude_south, translate_origin, 
				    camtype));
  }
  else if (p.defined("camera_position")) {
    // only one camera defined free-form
    vector<vec3d> positions, directions, up;
    vector<T_float> fovs;
    positions.push_back(p.getValue("camera_position", vec3d ()) );
    directions.push_back(p.getValue("camera_direction", vec3d ()) );
    up.push_back(p.getValue("camera_up", vec3d ()) );
    const T_float fov = p.getValue("camerafov", double ());
    fovs.push_back(fov);
    e.reset(new full_sed_emergence (positions, directions, up,
				    fovs, npixels, translate_origin, camtype));
  }
  else if (p.defined("camera_positions")) {
    // Read camera positions from file. This can have two different
    // formats: either theta and phi on a line, or complete 3-vectors
    // of position, direction, and up, and in addition the fov.
    const string file_name = 
      word_expand (p.getValue("camera_positions", string ())) [0];
    ifstream file (file_name.c_str());
    if (!file) {
      cerr << "Error: Can't open camera positions file "
	   << file_name << endl;
      exit (1);
    }
    else {
      cout << "Reading camera positions from file " << file_name << endl;
    }

    int npos=0;    
    {
      // Figure out the file format (this would be SO easy in Python...)
      ifstream tmpfile(file_name.c_str());
      string line;
      getline(tmpfile, line);
      istringstream is(line);
      T_float d;
      while(is >> d) {
	//is >> d;
	++npos;
      }
      cout << npos << " columns in camera positions file" << endl;
    }

    if(npos==2) {
      // theta and phi
#ifdef BROKEN_ISTREAM_ITERATOR
      // this statement is too fragile...
      istream_iterator<float_pair> a(file);
      istream_iterator<float_pair> b;
      std::vector<float_pair> positions(a,b);
#else
      std::vector<float_pair> positions(istream_iterator<float_pair>(file),
					(istream_iterator<float_pair>()));
#endif
      cout << positions.size() << " camera positions defined." << endl;

      // Convert to radians
      std::transform (positions.begin(), positions.end(), positions.begin(),
		      bind2nd(multiplies<float_pair>(),
			      float_pair(constants::pi/180,constants::pi/180)));
      const T_float distance = p.getValue("cameradistance", double ());
      const T_float fov = p.getValue("camerafov", double ());
      e.reset(new full_sed_emergence (positions, distance, fov, 
				      npixels, translate_origin, camtype ));
    }
    else if (npos==10) {
      // pos, dir, up, fov
      vector<vec3d> pos,dir, up;
      vector<T_float> fov;
      while(!file.eof()) {
	vec3d p,d,u;
	T_float f;
	file >> p[0] >> p[1] >> p[2]
	     >> d[0] >> d[1] >> d[2]
	     >> u[0] >> u[1] >> u[2]
	     >> f;
	if(!file.eof()) {
	  pos.push_back(p);
	  dir.push_back(d);
	  up.push_back(u);
	  fov.push_back(f);
	}
      }
      cout << pos.size() << " camera positions defined." << endl;
      e.reset(new full_sed_emergence (pos, dir, up, fov, npixels, 
				      translate_origin, camtype ));
    }
    else {
      cerr << "Error: Can't understand what the camera positions file contains." << endl;
      exit(1);
    }
  }
  else {
    cerr << "Error: Can't understand how you've defined cameras!" << endl;
    exit (1);
  }

  e->write_parameters(file, units);
  open_HDU (file, "MCRX").addKey("N_CAMERA",e->end()- e->begin(),
				 "Number of viewpoints");
}


void mcrx::Mcrx::setup_units () 
{
  const string output_file_name = 
    word_expand (p.getValue("output_file", string ())) [0];

  auto_ptr<FITS> output;
  output.reset(new FITS (output_file_name, Read));

  ExtHDU& pd_hdu = open_HDU(*output, "MCRX");
  pd_hdu.readKey("lengthunit", units["length"]);
  pd_hdu.readKey("timeunit", units["time"]);
  pd_hdu.readKey("massunit", units["mass"]);
}

/** Writes a history blurb to the specified FITS HDU.  This contains
    information such as final version, time, user, host, etc. */
void mcrx::Mcrx::write_history (HDU & output) const
{
  ostringstream history;
  history << "This file was created by Sunrise (mcrx) version "
	  << mcrx::version() << ", on ";
  time_t now = time (0) ;
  char*now_string = ctime (& now);
  now_string [24] = 0;
  history << now_string << " by user ";
  const char* u=getenv ("USER");
  assert(u != 0);
  string user (u);
  const char* h=getenv ("HOST");
  if(h == 0)
    h = getenv ("HOSTNAME");
  if(h == 0) {
    cerr << "Can't determine host name, HOST or HOSTNAME not set.\n";
    h="UNKNOWN";
  }
  string host (h);
  history << user << " on host "<< host << ".";
  //cout << history.str();
  output.writeHistory(history.str());
}


/// Sleeps, looking for termination file once per minute.
void mcrx::Mcrx::file_sentry::operator () ()
{
  const string termination_file ("mcrx_terminate");
  
  cout << "File sentry started" << endl;
  int fd;
  while (!self->terminator()) {
    //cout << "file sentry trying\n";
    if ((fd = open (termination_file.c_str(),O_RDONLY, mode_t ())) >= 0){
      // file found -- KILL!
      close (fd);
      //unlink (termination_file.c_str());
      cout << "File sentry found termination file -- terminating!\n";
      self->kill();
      return;
    }

    // sleep
    struct timespec sleepy_time;
    sleepy_time.tv_sec = 60;
    sleepy_time.tv_nsec = 0;
    nanosleep (&sleepy_time, 0);
  }
}


/** Constructor creates a time_sentry with the specified margins and
    connected to the specified Mcrx object.  */
mcrx::Mcrx::time_sentry::time_sentry (Mcrx* m, int cm, int wm):
  self (m), start (time (0)), cpu_margin (cm), 
  wall_margin (wm), dump_on_exit(true)
{
  const char*l = getenv("WALL_CLOCK_LIMIT");//
  if (l){
    std::istringstream ist (l);
    ist >> wall_limit;
  }
  else{
    wall_limit = 0;
    cout << "wall clock limit not found, check disabled" << endl;
  }
}

/// Sleeps, checking wall clock and CPU time limits once per minute.
void mcrx::Mcrx::time_sentry::operator () ()
{
  struct rusage usage;
  struct rlimit limit;
  getrlimit (RLIMIT_CPU, &limit);
  int cpu_limit =limit.rlim_cur;
  cout << "Time sentry started, CPU limit:" << cpu_limit << "s, wall clock limit " << wall_limit<< "s" << endl;
  
  while (!self->terminator()) {
    getrusage (RUSAGE_SELF, &usage);
    getrlimit (RLIMIT_CPU, &limit);
    int cpu_usage =usage.ru_utime.tv_sec+usage.ru_stime.tv_sec;
    int wall_usage = time (0) - start;

    if ((cpu_limit > 0) && (cpu_limit - cpu_usage < cpu_margin)) {
      // we are running out of CPU time
      cout << "Time sentry:running out of CPU time -- restarting!\n";
      self->restart();
      return;
    }

    if ((wall_limit > 0) && ( wall_limit - wall_usage < wall_margin)) {
      // we are running out of wall clock time
      cout << "Time sentry:running out of wall clock time -- restarting!\n";
      self->restart ();
      return;
    }
      
    cout  << "CPU left: "<< cpu_limit - cpu_usage <<  "\tWall clock left " << wall_limit - wall_usage << endl;
    
    // sleep
    struct timespec sleepy_time;
    sleepy_time.tv_sec = 60;
    sleepy_time.tv_nsec = 0;
    nanosleep (&sleepy_time, 0);
  }
}

/** Prints used CPU and wall clock time. */
void mcrx::Mcrx::time_sentry::print_used() const
{
  struct rusage usage;
  getrusage (RUSAGE_SELF, &usage);
  int cpu_usage =usage.ru_utime.tv_sec;
  int cpu_susage =usage.ru_stime.tv_sec;
  int wall_usage = time (0) - start;

  cout << "\nUsed resources:"
       << "\n\tCPU time: usr " << cpu_usage
       << " s\n\t          sys " << cpu_susage
       << " s\n\t          tot " << cpu_susage+cpu_usage
       << " s\n\tWall time: " << wall_usage << " s" << endl;
}

    
/** The main function for the mcrx executable.  The functionality is
    in the Mcrx object, all this function does is figure out which
    function to call and when.  It also sets up the sentries.
    \ingroup mcrx */
int main (int argc, char** argv )
{
  using namespace mcrx;

  cout << "mcrx (Sunrise) version " << mcrx::version() << ". Copyright 2004-2011 Patrik Jonsson.\nThis program comes with ABSOLUTELY NO WARRANTY. Sunrise is free software, and\nyou are welcome to redistribute it under certain conditions. See the GNU\nGeneral Public License for details.\n\n";

  if (argc != 2) {
    cout << "Usage: mcrx <config file>\n";
    exit (1);
  }

  assert(blitz::isThreadsafe());
  Preferences p;
  p.readfile (argv [1]);

  p.setValue("pwd", string (getenv ("PWD")), "Working directory") ;
  p.setValue("configfile", string (argv [1]) , "Configuration file") ;

  // check if hpm should be enabled (default is disabled)
  if (p.getValue("use_hpm", bool (false), true)) {
    hpm::init();
    hpm::enable ();
  }
  else
    hpm::disable ();

  // check if counters should be disabled (default is enabled)
  if (p.getValue("use_counters", bool (true), true))
    counter::enable = true;
  else
    counter::enable = false;

  // check if CCfits should be verbose (default is not)
  CCfits::FITS::setVerboseMode (p.getValue("CCfits_verbose", 
					   bool (false), true));

  numa_info();
  if(p.getValue("interleave_common_data", false, true))
    interleave_allocations();

#ifdef WITH_CUDA
  // set up CUDA if it's used
  if(!mcrxcuda::cuda_init(p.defined("cuda_device") ? 
			  p.getValue("cuda_device", int()) : 0)) {
    // CUDA init failed, turn it off
    p.setValue("use_cuda", false);
  }
#endif

#ifdef WITH_AREPO
  // do necessary initialization for Arepo
  arepo::init(&argc, &argv);
#endif

  if (p.defined("mcrx_data_dir"))
    particle_base::inp_dir = 
      word_expand (p.getValue("mcrx_data_dir", string ())) [0];
    
  Mcrx theRun (p);

  Mcrx::file_sentry fs (&theRun);
  Mcrx::time_sentry ts (&theRun, p.getValue("cpu_time_margin", int ()),
			       p.getValue("wall_clock_margin", int ()) );
  boost::thread_group sentries;
  sentries.create_thread(fs );
  sentries.create_thread(ts );
  
  if (!theRun.terminator())
    theRun.initialize_output_file();

  theRun.setup_units();

  theRun.execute_stages();

  // first check: if skip_postprocessing is set, then end right here,
  // right now, with a successful exit code
  if (!theRun.terminator() &&
      p.defined( "skip_postprocessing") &&
      p.getValue("skip_postprocessing", bool ()) ) {
    cout << "Skipping postprocessing, we're done." << endl;
    //hpm::hpmterminate (0);
    return 42;
  }
  
  // if we have defined a separate submit command for postprocessing,
  // we end here and submit that.
  if (!theRun.terminator()&&p.defined( "postprocess_command_variable")) {
    const char*rsc = getenv (p.getValue("postprocess_command_variable", 
					string()).c_str());
    if (rsc) {
      string resubmit (rsc);
      // execute resubmit command
      cout << "Executing postprocess command: " << resubmit << endl;
      system (("exec " + resubmit).c_str());
      //hpm::hpmterminate (0);
      return 42;
    }
    else {
      cout << "No command defined for postprocessing, proceeding with postprocessing now." << endl;
    }
  }

  // postprocessing stages
  theRun.calculate_attenuated_images();
  theRun.generate_integrated_SED();
  //theRun.add_infrared_templates();
      
  //hpm::hpmterminate (0);

  // should we resubmit
  if ((theRun.terminator() == terminator::restart_) &&
      p.defined("resubmit_command_variable")) {
    const char*rsc = getenv (p.getValue("resubmit_command_variable", 
					string()).c_str());
    if (rsc) {
      string resubmit (rsc);
      // execute resubmit command
      cout << "Executing resubmit command: " << resubmit << endl;
      system (("exec " + resubmit).c_str());
      return 66;
    }
    else {
      cout << "No command defined for resubmission" << endl;
      return 102;
    }
  }

  // set appropriate exit value
  if (theRun.terminator())
    return 7;

  cout << "mcrx complete.\n";

  return 0;
}


/** This function selects between the Arepo and octree code paths. */
void mcrx::Mcrx::execute_stages()
{
  const string output_file_name = 
    word_expand (p.getValue("output_file", string ())) [0];
  auto_ptr<FITS> output (new FITS(output_file_name, Read));
  string nbodycod;
  ExtHDU& sh = open_HDU(*output, "SFRHIST");
  sh.readKey("NBODYCOD", nbodycod);

  if (nbodycod=="arepo") {
#ifdef WITH_AREPO
    cout << "Selecting Voronoi mesh code path\n";
    string snapshot_file;
    string parfile;
    sh.readKey("snapshot_file", snapshot_file);
    snapshot_file = strip_suffix(snapshot_file);
    sh.readKey("simparfile", parfile);
    const int n_threads = p.defined("n_threads")?
      p.getValue("n_threads", int ()): 1;
    arepo::load_snapshot(snapshot_file,parfile, n_threads, units);
    //MUST close the file here otherwise we can't open it ReadWrite in the stages
    output.reset();
    execute_arepo_stages();
#else
    cerr << "Not compiled with Arepo support" << endl;
    exit(1);
#endif
  }
  else {
    cout << "Selecting octree mesh code path\n";
    output.reset();
    execute_adaptive_stages();
  }
}
